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ABSTRACT 

We present an ab initio approach to the solar coronal heating problem by modelling a small part 
of the solar corona in a computational box using a 3D MHD code including realistic physics. The 
observed solar granular velocity pattern and its amplitude and vorticity power spectra, as reproduced 
by a weighted Voronoi tessellation method, are used as a boundary condition that generates a Poynt- 
ing flux in the presence of a magnetic field. The initial magnetic field is a potential extrapolation 
of a SOHO/MDI high resolution magnetogram, and a standard stratified atmosphere is used as a 
thermal initial condition. Except for the chromosphcric temperature structure, which is kept fixed, 
the initial conditions are quickly forgotten because the included Spitzer conductivity and radiative 
cooling function have typical timescales much shorter than the time span of the simulation. After 
a short initial start up period, the magnetic field is able to dissipate 3 — 4 x I0 6 ergs cm" 2 s _1 in a 
highly intermittent corona, maintaining an average temperature of ~ 10 6 K, at coronal density values 
for which emulated images of the Transition Region And Coronal Explorer(TRACE) 171 and 195 A 
pass bands reproduce observed photon count rates. 

Subject headings: Sun: corona - Sun: magnetic fields - MHD 
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1. INTRODUCTION 

The heating mechanism at work in the solar corona 
has puzzled researchers for more than six decades. Sev- 
eral heating mechanisms have been proposed. One 
of the earliest models involved accretion of interstel- 
lar matter coupled with convection in the chromosphere 
( Bond i"e"t al.l ^947) , but with improved measurements 
two main groups of heating models emerged: Wave 
(AC) heating and electric current (DC) heating. These 
mechanisms move energy from the photospheric kinetic 
energy reservoir to the internal heat reservoir in the 
corona through th e magnetic fi eld. AC heating (al- 
ready proposed by lAlfvenl 11947ft depends on the mag- 
netic field being moved around in the solar photosphere 
faster than the disturbances can propagate through the 
whole magnetic loop, i.e. faster than the Alfven cross- 
ing time. Only torsional Alfven waves can reach the 
corona, while magneto-sonic wave modes are diffracted 
and dissipated due to the strong wave speed gradient in 
the chromosphere/transition region. Waves have been 
detected both in the near Sun corona (se e for instance 
IAschwanden1ll987t lOfman fc Davilalll997j) as well as in 
the solar wind b y a number of spacecrafts su ch as the 
Helios satellites l)Neubauer fc Musmann 1977J) so there 
is firm evidence for their existence, but not of their 
ability to supply enough heat. Alfven waves do not 
easily d issipate in the corona in the absence of phase 
mixing ( Hevvaerts fc Priest! 1983(1 or resonant absorption 
i|Davilalll987|) . which require strong phase speed gradi- 



ents to be efficient. 

DC heating relies on motions in the photosphere 
that change on time scales longer than Alfven cross- 
ing times, making the magnetic field remain close to 
an equilibrium state. The DC heating mechanism 
dissipates ene rgy through conventional Joule heating 
(proposed by Parker} 11973). or rec onnection heating 
(proposed by iClencross et alJ 11974 based on X-ray 
observations) dissipated through flares of all sizes, 



including nano-flares, or through a hierarchy of current 
sheets that may encompass bo th of the above processes 
(|Galsgaard fc Nordlundl fl996). Until recently, nano- 
flares appeared to be the most promising candidate, 
but it now seems that the energy dissipated in the 
observed flare distribution, extrapo lated to a cut-off 
at nan o-flare energies, is too s mall l|As chwanc}en e^ alJ 
2jM IParnell fc .TirnTI 1200?! lAsrhwauden fc Pa.rnelll 
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20021 lAschwanden fc Charbonneaul 12002(1 DC heat- 
ing has received considerable attention since i t was 
propo se d some thirty yea rs ago (e. g. iParkerl 119721 
1983} iSturrock fc Uchidal Il981t Ivan Ballegooiier 
19861 iMikic et alJ Il98fl iHevvaerts fc Priestl I199S 



Lons-cone fc Sudar 



S 



Galsgaard & NordlunJ E 



Hendrix et al]ll99fitlGomez et al.ll200(l to mention just 
a few). It has been established that the DC heating 
mechanism is feasible, but not to what extent it can 
provide enough energy to heat the solar corona. 

The main problem in modelling solar dissipative pro- 
cesses is the length scales involved. The dissipative 
length scale is of the order of 1 m for e.g. a flare 
with a time scale of the order 100 s, while the large 
scale magnetic field distribution has a scale of many 
Mm, making fully resolved 3D simulations of such pro- 
cesses pra ctically i mpos sible. However, the theoretical 
results of IParkerl fl!97S ) and the numerical results o f 
IGalsgaard fc Nordlundl l(1996|) and lHendrix et all (|1996Tl 
have shown that it is not necessary to resolve the dis- 
sipative length scales because the total dissipation only 
depends very weakly on the resolution, and if anything 
increases with increasing resolution. Results from low 
resolution modelling of dissipative processes may thus 
be used to obtain estimates, or at least establish lower 
limits on the dissipated energy. This provides the ba- 
sis for trying to create numerical ab initio simulations of 
the solar corona, with essentially no free parameters. As 
discussed in more detail below, the initial and boundary 
conditions that we use are minimal assumptions, in the 
sense that including more detailed conditions would pro- 
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Fig. 1. — Radiative cooling function (full) and implemented cool- 
ing function (dashed) with quenched cooling in the photospheric 
and chromospheric temperature regimes 



duce more heating. This allows us to give a firm answer 
to whether DC heating can provide the energy required 
and thus which heating mechanism is the principal one 
in the solar corona. 

Initial results of t his e ffort were presented in 
Gudi ksen fe Nordlundl l)2002(l and here we present fur- 
ther developments including a realistic photospheric ve- 
locity driver and a non uniform mesh, making it possible 
to resolve the density drop in the transition region and 
expand the simulated volume. 

2. MODEL 

We solve the fully compressible MHD equations on a 
non-uniform staggered mesh, in the form 



dp 
dt 
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dt = -^- eu 



PV -u 

— V • Pspitzcr + Qvisc + Qjouic — n ion n e A , (6) 

where p is the mass density, u the velocity vector, r the 
viscous stress tensor, P the gas pressure, J the electric 
current density, B the magnetic flux density, g the grav- 
itational acceleration, p the vacuum permeability, E the 
electric field strength, ij the magnetic diffusivity, e the 
internal energy per unit volume, P snitzer is the energy 
flux due to the Spitzer conductivity l)Spitzerlll956j) along 
the magnetic field, Qvisc is the viscous heating, Qjouio is 
the Joule heating, and A is the cooling function for the 
optically thin coronal plasma with rii on and n e being the 
number density of ions and electrons. 

The code is explicit and, because of the staggered mesh 
method, conserves the magnetic divergence to machine 
precision. The staggered mesh method means that the 



physical variables are not aligned in space. This is an ad- 
vantage for the derivative operators, but has the down- 
side that when physical variables at different locations 
are involved interpolation is needed to realign the vari- 
ables. The code uses 6'th order derivative operators and 
5'th order interpolation operators, together with a 3'rd 
order Runge-Kutta time stepping method with variable 
time step. The code has been extensively subjected to 
several standard tests as well as used in several physi- 
cal regimes, for instance 3-D tur bulence and magneto- 
convection llNordlund et aflll994D . magnetic dissipation 
llGalsgaa nij^Nordlund jj.99^1. buoyant magnetic flux 
tubes ( Dorch & Nordlund 1998), em ergence of flux rope s 
through the solar convection zone llDorch et alj[l 999). 
and supersonic MHD-turbulence l|Padoan fe Nordlundl 
1999). 

The computational box of 60 x 60 x 37 Mm 3 is resolved 
on a grid of 150 3 grid points where the vertical scale is 
non-uniform. The vertical axis starts at the photosphere, 
with high resolution (~ 0.15 Mm) in order to resolve the 
small pressure scale height in the photosphere and chro- 
mosphere. Above the transition region the vertical reso- 
lution gradually approaches 0.25 Mm. 

2.1. Radiative cooling 

The cooling function includes the most important tran- 
sitions for H, He, C, O, Ne, Fe and bremsstrahlung, 
with a cut off at low temperatures to avoid catas- 
trophic cooling in the chromosphere (see Fig. |y. The 
cooling function is based on the ionisation and recom - 
bination rates given bv lArnaud fc Rothenflud l)1985|k 
IShull fc van Steenberel J1982D and using the collisional 
excitation rate s found through the HAO -DIAPER Atom 
data package ij.Tudge fc Meisnerl 119941) . The ions are 
treated by assuming ionisation equilibrium making it 
possible to derive radiative losses as a function of electron 
temperature. 

Further down in the atmosphere, where the radiation is 
optically thick for at least some wavelengths, the cooling 
processes cannot be treated with a simple cooling func- 
tion. The cooling function shown in Fig. ^full) uses a 
cutoff proportional to exp (— r) with the optical depth 
r oc Pgas- In this simulation, the lower part of the 
chromosphere and the photosphere are instead kept near 
a fixed average temperature profile using a Newtonian 
cooling mechanism that forces the local temperature to- 
wards a preset function of height on a timescale of about 
0.1 s in the photosphere. The time scale increases with 
decreasing density as p~ 167 , and the effect thus becomes 
negligible in the corona. 

2.2. Boundary conditions 

The vertical boundary conditions arc handled by using 
ghost zones, while the box is periodic in the horizontal 
directions. The density is extrapolated into the ghost 
zones, and to prevent heat from leaking into or out of 
the box, the temperature gradient is set to zero at the 
upper boundary. The upper boundary is relaxed towards 
the average temperature at the upper boundary with a 
typical time scale of about 30 s. The vertical velocity 
is zero on the boundaries, while the horizontal velocity 
on the lower boundary is set by a driver (See ^2.3fl . At 
the upper boundary the vertical gradient of the horizon- 
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Fig. 2. — Velocity amplitude versus spherical harmonic 
wavenumber for the three layer granulation model (full) with the 
error bars showing maximum deviation through a period of 1 hour. 
The dashed line shows u oc k. 



tal velocity is set to zero and we use a potential field 
extrapolation for the magnetic field in the ghost zones. 

2.3. The Voronoi tessellation driver 

The solar photospheric velocity field is characterised by 
granular structures with a range of sizes, from granules 
with a typical size of ~ 1 Mm to super -granules of size 
~ 20—30 Mm and even bigg er giant cells l)Simon fc Weissl 
119681 lHathawav et al.l l2000'). These structures have typ- 
ical velocities that scale with wavenumber roughly as 
u cx k and have turnover times of the order r cx fc~ 2 , 
all the way from beyond super-gran ular scales to the 
velo ci ty maximum at granu l ar scales iStein fc Nordlundl 
Il998t IHathawav et alJl200Cl IShine et alJl2000D . The ve- 
locity field injects energy into the magnetic field through 
a Poynting flux, making a realistic reproduction impor- 
tant. 

Our goal was to create a driver that would match these 
statistical features of the solar velocity spectrum, as well 
as the geometrical features of the granular pattern. The 



velocity driver is partially bas ed on the theory of Voronoi 
tessellation (see for instance lOkabe et al.lll992T) specifi- 
cally mul tiplicatively weighted V oronoi tessellation, in- 
spired bv lSchriiver et al.1 l[r997a1 who showed that this 
method gives a good fit to observations of the granula- 
tion pattern. The theory applied to this case provides a 
way to split a 2-D surface into tiles, each representing a 
granule. To do this one has to chose a number of gen- 
erator points, each with a weight Wi (where i represents 
a specific tile) positioned at a;,-. Tile i then occupies an 
area where the following inequality is satisfied, 



> 



W j 



for all j^i 



(7) 



To include the intergranular lanes, we increased the 
complexity of Eq[7| 
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When both conditions are true x is in an intergranu- 
lar lane, and where only the left part of the inequality 
is true x is in granule i. h g is a number between zero 



and one and can be adjusted to mak e the intergranular 
area roughly 33-45% of the total area l|Stein fc Nordlundl 
1998) based on velocity maps. The relative area of the in- 
tergranular lanes can be estimated by assuming perfectly 
round granules, 
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2.3.1. The initial tessellation 



If the generator points were placed randomly there 
would be both unusually large and unusually small gran- 
ules generated, because the density of generator points 
would be low in some places and high in others. Since ob- 
served granules are rather uniform in size, only a certain 
range of sizes of fully grown granules can be accepted. 
This places restrictions on the placement of the genera- 
tor points. The tessellation is built up by keeping a list 
of the weights at each position, independent of whether 
it is in an intergranular lane or not, updating the list as 
more granules are placed. In this way a weight function 
with a value at each position is generated, 



W(x) = ■ 



Wi 



> 



where 

W 3 



for all j ^ i . (10) 



To generate a new granule, its generator point is placed 
where the weight function has a value lower than a 
certain threshold Wthrcsh- This allows more generator 
points to be placed only until the weight function is above 
the threshold at all positions. By choosing an appropri- 
ate threshold carefully the density of generator points is 
such that granules of the right size are created. 

The velocity vector points radially away from the gen- 
erator points inside the granules, and along the gran- 
ular boundaries in the intergranular lanes. The hor- 
izontal velocity profile i nside meso-granules was mod- 
elled bv lSimon fc Weiss (1989) by fitting a simple model 
to observations done with the SOUP filter on Spacelab 
2. The horizontal velocity profiles were successfully fit- 
ted with a radial velocity profile of the form v*h(r) = 
r exp (— (r/r n ) 2 ). We found a better fit to the solar veloc- 
ity power spectrum by using Vh(r) — r 2 exp (— (r/rrj) 2 ), 
giving a flatter profile in the centre of the granules. In 
the intergranular lanes the velocity is a sum of the veloc- 
ity of the two granules adjacent to the intergranular lane. 
Initially this creates stagnation points in the intergran- 
ular lanes, on lines connecting the generator points, but 
such stagnation points disappear when multiple velocity 
fields are added (see H2. 3.2(1 . 

2.3.2. Time evolution 

The tessellation is made time dependent by making the 
weights of the generator points time dependent. We have 
chosen to give the weights of the generator points a time 
dependence specified by 



Wi (t) = Wiexp - 



t - 1 



0,i 



n 



(ii) 



where Wi is a constant giving the maximum weight, tg^ 
is the time of maximum weight, and q is an even integer, 
indirectly controlling the growth time of the granule. We 





a id 



■■ I ■■ — . 1 1 : n . .i? - . — r. . .' 

2D Sfl 40 



so kg 



jo iD 

X [ Mm ] 



fib 6u 



"0 



2D 30 40 
H [ Mm ] 



Fig. 3. — Three snapshots of 3000 corks released in the simulated velocity field, without magnetic velocity quenching, at t 
(left), t = 15 min (middle), and t = 45 min (right) 
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bring the granule to life when Wi j dx > Wthresh where dx 
is the grid resolution, which effectively sets to^ when 
is chosen. When t > t oi the granule is decaying and it 
is removed when Wi/dx < Wthresh- 

When the weight function W(x) falls below the thresh- 
old in a large connected region, each new granule covers 
only a small area, causing many small granules to ap- 
pear, thus making the area crowded with small gran- 
ules of approximately equal size. This problem is solved 
by creating a binary function A{x, t) , which is set to 
one within 80 % of the expected radius for granule i at 
t = to,i as long as t < to.i and otherwise zero. This re- 
gion of avoidance keeps new granules from being formed 
in the immediate vicinity of other newly formed gran- 
ules. Therefore there are two criteria a new generator 
point (n) must fulfil when it is generated, in order to be 
part of the tessellation, 

res ii 

A(x n ,t)=0 . 

This ensures that there is room for a new granule and 
avoids overpopulation of new granules later when the 
granules are growing. At the same time the weight func- 
tion evolves such that it makes new granules appear only 
in intergranular lanes, because this is the location where 
the weight function first drops below the threshold. 

In order to create the velocity-scale relation of the Sun, 
several granular patterns with different scales must be 
superposed. This involves creating a number of indepen- 
dent granulation patterns with different typical scales as 
described above and adding the velocities. We have cho- 
sen three layers, with typical granular radii of 1.3, 2.5, 
and 5.1 Mm respectively, with velocities and lifetimes 
following the observed scaling relations. This choice cor- 
responds to sizes ranging from those of large granules to 
the size of mesogranules, and produces on average 880, 
280 and 80 granules respectively, in the three size groups. 
The smallest scale is set by the resolution and the largest 
scale is set by the box size used in the simulation. 

Creating the velocity field in this way makes the ve- 
locity field ordered on many scales, but also makes it 
too laminar. The vorticity spectrum from the Voronoi 
driver is too we a k com pared to the simulations of 
iStein fc Nordlundl l)1998j) . This has been corrected by 



using a Helmholtz projection to separate the rotational 
and irrotational part of the velocity field and then am- 
plifying the ro tational part to the sam e level as in the 
simulations bv IStein fc Nordlundl l)1998ft . 
The velocity is quenched by a factor 

_ i + r 2 , 12 , 

J quench — ^ ^ 3|3~ 2 ' 

depending on the plasma (3 — P ga s/ Phiag- In the 
strongest magnetic field regions the quenching reduces 
the velocity amplitude by ~ 60 %, in order to reproduce 
the magnetic fields ability to quench convective motions 
in the Sun. 

The velocity - scale relation of the driver, averaged 
over 45 min, is shown in Fig. |5J To illustrate the time 
development of the driver, snapshots of corks flowing in 
the time-dependent velocity field produce by the driver 
are shown in Fig. [21 

2.4. Initial conditions 

Two initial stratifications w ere used. The first was 
taken from the FAL-C model ijFontenla et alJ IT993fl in 
the photosphere, chromosphere and transition region, 
while in the corona a simple isothermal profile with 
T = 10 6 K was use d (see FiglH. The second was in- 
spired by the work of ICarlsson fc Steinl (see for instance 
ICarlsson fc Stemll2002F who argue for a highly time de- 
pendent, low mean temperature chromosphere with no 
average temperature incr ease. This model has been crit- 
icised by Kalkofen ( 2001), who argues for a warmer chro- 
mosphere. On the other hand observations of CO by 
lAvresI l)2002l) point toward s an even cooler ch romosphere 
than the one proposed bv ICarlsson fc Steml . The choice 
of chromospheric model turns out to have only a minor 
effect on the coronal dynamics and the coronal heating. 

To keep the number of free parameters at a minimum 
we chose to make the initial magnetic field a potential 
extrapolation of an observed active region on the Sun. 
We chose to scale down a high-resolution SOHO/MDI 
observation of active region 9114 from August 8, 2000. 
The cropped data spans 500 x 500 pixels corresponding 
to roughly 225 x 225 Mm and was therefore scaled down 
to fit in the computational box of 60 x 60 Mm. The dis- 
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Fig. 4. — Initial stratification for density (full), gas pressure 
(dashed) and temperature (dotted). Also shown at the top is the 
vertical-scale, with a small line at each grid point 
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Fig. 5. — Horizontally averaged energy deposition from resistive 
dissipation (full), Spitzer conductivity (long dash), convective flux 
(dash dotted) and radiative cooling (dotted), as a function of height 



tribution of flux on the solar surface seems to be approx- 
imatel y self-similar down to scales ^180 km in the net- 
work (Schriiver et alJll997bl iWang et al.lll995j) . smaller 
than our resolution. If this also holds in active regions, as 
seems to be the case ijMeunierl 12003). this would justify 
our "rescaling" of the MDI magnetogram. It is however 
necessary to point out that our computational box does 
not hold an average size solar active region, but at most 
a small active region, and therefore direct comparisons 
of these results with larger size active regions should be 
done with care. 

3. RESULTS 

The simulation goes through a roughly five minute 
start-up period, where neither thermal conduction nor 
radiative cooling are turned on. This is related to the 
magnetic field initial condition. Since the initial mag- 
netic field is potential it is in a minimum energy state 
and is therefore initially not able to dissipate any en- 
ergy. In the first few minutes, thermal conduction and 
radiative cooling would therefore cool the corona mono- 
tonically, and make the atmosphere tend to collapse. 
After the start-up period, when the magnetic field has 
reached a dissipative state, the radiative cooling and the 
thermal conduction is turned on. After another approx- 
imately five minutes solar time a statistical equilibrium 
is reached, with a nearly time independent temperature 
distribution. 

3.1. Energy Balance 

The thermodynamic equilibrium of the atmosphere is 
difficult to quantify since the energy fluxes are very inter- 
mittent, but several conclusions can still be drawn from 
a horizontal average of the energy flux divergences, i.e. 
the deposited energy. The energy balance is shown in 
Fig. [5] This balance is for a typical point in time, and 
it is obvious that the heating in the corona is decreasing 
with heig ht, as found from TRACE data by a number of 
autho rs ijSchriiver et al.lll999t lAschwanden et alJ feOOOa. 
1200 1|) . The detailed height dependence of the heating 
in individual loops is invest igated in a separate paper 
llCxudiksen fc Nordlundll200l . 

It is noteworthy and interesting that the convective 



flux redistributes as much energy as the Spitzer conduc- 
tivity does just above the transition region. The heat- 
ing provided by the magnetic dissipation is transported 
by the Spitzer conductivity downwards from the transi- 
tion region to the upper chromosphere, where the Spitzer 
conductivity is no longer effective due to the low temper- 
ature. The convective flux removes energy at the same 
location because of evaporation of chromospheric mate- 
rial into the corona. The energy loss in the lower corona 
and transition zone due to the Spitzer conductivity is 
balanced by the heating, which dominates the Spitzer 
conductivity at larger heights. The optically thin radia- 
tive losses go to zero in the lower chromosphere and pho- 
tosphere due to the quenched radiative cooling function 
that we adopt. The heating goes to very high values, but 
in spite of that the temperature stays at chromospheric 
temperatures due to the Newtonian cooling with which 
we represent the complex chromospheric radiation losses 
at these height. It is not possible to estimate whether any 
further heating is needed in the chromosphere until we 
can consistently treat both convective and radiative en- 
ergy fluxes through the lower boundary, as well as radia- 
tive losses in the optically thick lower atmosphere. There 
should be no major effect on the conclusions drawn about 
the heating of the transition region and corona due to this 
inadequacy, because we hold the temperature structure 
in this region of the Sun close to what it is deduced to 
be observationally. Incorporating both enthalpy flux and 
an appropriate treatment of radiation in a near ab-initio 
chromospheric model will take a considerable effort, and 
is a challenge for the future. 

3.2. Magnetic field configuration and the dissipative 
heating 

The energy supplied to the corona originates from the 
Poynting flux through the lower boundary. In general, 
magnetic dissipative heating is generated from stressing 
of the magnetic field. The stressing is largest in the 
photosphere and low chromosphere, where the plasma 
/3 is larger than or of the order of unity in a significant 
fraction of the volume. At these low heights the mag- 
netic field still has a very intermittent structure, reflect- 
ing the intermittent distribution set up by fluid motions 
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in the photosphere and below (represented by the ini- 
tial and boundary conditions in our numerical model). 
The plasma (3 is low in regions with high magnetic field 
strengths, but is still high in intervening regions with low 
magnetic field strengths. 

In the high chromosphere and above the magnetic field 
is nearly space-filling and the plasma (3 is less than unity 
almost everywhere (with the exception of small neigh- 
bourhoods around magnetic null points). The magnetic 
field thus dominates the dynamics and tends to be 'nearly 
force-free', in the sense that the Lorentz force is small and 
the electric currents tend to be nearly parallel with the 
magnetic field, 

J « aB. (13) 

As is well known and easily proven, a must be constant 
along each magnetic field line in a truly force- free field, 
but can certainly vary from one field line to another. 
Also, since a dynamically evolving coronal magnetic field 
cannot be entirely force-free, a should not be expected 
to be exactly constant even along magnetic field lines. 

In order to illustrate the variability of a over horizontal 
planes, Fig.HJshows images of the quantity | JB | w aB 2 
(we prefer to show this quantity rather than | J ■ B/B 2 |, 
which would tend to be dominated by regions with low 
B). 

If a was constant these images would effectively be the 
same as images of B 2 . The actual images show that a 
varies significantly over horizontal planes, that it varies 
most rapidly at low heights, and that it certainly cannot 
be approximated with a constant at any height. The 
larger horizontal scale of the variability at larger heights 
is due to the fact that field lines (along which a is at least 
approximately constant) fan out as they reach higher, 
and hence spread the variability of a over larger patches. 

We next comment on the absolute magnitude of a, and 
its consequences for the overall field distribution. Dimen- 
sionally, a is an inverse length; the length along a field 
line over whic h it twists around itse l f. As previously 
established by Galsgaard & Nordlund (1996), field lines 
that are stressed by a boundary typically are twisted only 
about once from end to end, corresponding to values of 
a of the order of l/L, where L is the field line length. 
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Fig. 7. — PDF of the electric current squared as a function 
of height, dark colours mean higher PDF with logarithmic colour 
table. Over-plotted is the horizontal average of the magnetic field 
squared (dotted) and the horizontal average of the current squared 
(dashed) 



As illustrated by Fig. El a is only correlated over small 
patches near the lower boundary. This implies that the 
perturbations of the magnetic field line directions must 
be small (except in special regions s uch as in the neigh- 
bourh ood of quasi-separatrix layers ijPriest &; Demoulhl 
1995)), and hence that the overall perturbations of the 
magnetic field strength distribution, relative to that of a 
potential field, must also be small. This means that to 
obtain an overall picture of the magnetic field strength 
distribution, a potential field extrapolation is sufficient 
in cases like this one, with no large scale shear. 

We expect the heating to generally be proportional to 
the electric current squared and if the field is in a non- 
linear force free state we also expect the electric current 
squared and thus the heating to be proportional to the 
magnetic field strength squared. Figure [7\ shows a con- 
tour plot of the distribution of heating with height in 
the atmosphere. One can see that the heating is largest 
at low heights due to the stressing of the magnetic field 
in a high j3 environment. Through the photosphere and 



7 



0.10 



o 

CL. 

'tn 



0.01 







■ 





















100 

Time scale [ s ] 

Fig. 8. — Dissipated energy as function of time scale in the 
upper chromosphere (full), transition region (long dashes), lower 
corona (dashes) and upper corona (dotted) in arbitrary units. The 
curves are normalised at the maximum time scale 
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Fig. 9. — Energy conversion efficiency as function of time during 
a 23 min period. The relative amount of energy converted into 
heating (solid) and the relative amount of energy lost to radiation 
(dashed) in the corona are shown 



chromosphere it declines by some four orders of magni- 
tude before reaching the transition region, above which 
the magnetic field is nearly force free and the heating 
is indeed approximately proportional to the magnetic 
field strength squared. Even though the electric cur- 
rent is highly intermittent in horizontal planes, the cur- 
rent relative to the magnetic field strength is almost con- 
stant along field lines above the transition region, as is 
to be expected when a is constant along magnetic field 
lines. Heating with th ese characteristics was proposed by 
iSchriiver et al.1 119991) and later deduc ed with data from 
TRACE bv lAschwanden et all (j2001^) . Further studies 
of heating models and observations ha ve all supported 
these heating function characteristics jM andrini et_aL| 
| 2000t iFolev et all I2QQ21 ISchriiver & Aschwandenl 120021 
Dcmoulin et alJl200.1 iSchmieder et al.ll2003L 

The time variability of the heating depends critically 
on height. Figure |S| shows on which time scales the en- 
ergy dissipation varies, calculated as k t \fP~ e where kt is 
the time-wavenumber, and P e is the power spectrum of 
the dissipated energy as a function of time at each point. 
This produces a plot showing on what timescales the en- 
ergy is dissipated. The energy dissipated is normalised 
to the total amount of dissipated energy. Without this 
normalisation the curves would be widely spaced along 
the y-axis, due to the exponential decrease in total dissi- 
pated energy with height. In general most of the en- 
ergy is dissipated as short time scale events. In the 
upper chromosphere, the energy comes primarily from 
events with a time scale shorter than the time resolu- 
tion of our snapshot data set (~ 38 s). Slightly further 
up into the corona, the fraction of energy dissipated at 
short time scales is less pronounced. The rapid events 
most likely correspond to "nano-flare" like reconnection 
events. With increasing spatial resolution, we would ex- 
pect to resolve smaller and smaller and hence increas- 
ingly more short-lived events. 

The energy supplied by the Poynting flux is distributed 
initially between changes of kinetic, thermal, potential 
and magnetic energy, and dissipation. After the system 
has settled into a statistical equilibrium the energy is 
only used to compensate for magnetic and kinetic dis- 



sipation. The energy is dissipated in the whole atmo- 
sphere, with the main part in the photosphere and chro- 
mosphere. The Poynting flux at the lower boundary is 
highly intermittent but the total Poynting flux is close 
to constant. Comparing the flux of energy through the 
lower boundary and the heating rate in the corona gives 
an overall energy conversion "efficiency" . The time evo- 
lution of the efficiency is shown in Fig. [^Jsolid) during a 
20 min period. During this time the numerical resistivity 
was changed several times, apparently without significant 
effect on the efficiency. The efficiency is in the range 5- 
10 %. That the changes in resistivity have little effect 
both on the Poynting flux into the system and on the ef- 
ficie ncy also confirms the results o f lHendrix et alJ l)199fif > 
and iGalseaard k, Nordluncfl l)1996|) . The radiative losses 
of the corona is in the range 1-3 % during the same time 
period (Fig.El dashed) making it less important than the 
Spitzer conductivity, as expected. 

The total amount of heating needed in active regions 
has been estimated several times based on UV and X- 
ray fluxes, and typically gives values in the range 10 6 - 
10 7 ergs cm" 2 s _1 . These values are in general hard to 
estimate precisely, since the total area of the active re- 
gion can be chosen based on different criteria. Figure 
ITU1 shows the average dissipated energy per second and 
square centimetre in our model, in gas that is emitting 
UV and X-rays (T > 5 x 10 5 K), in volumes below which 
the photospheric magnetic field is higher than a cer- 
tain threshold. If only regions with high magnetic field 
strength below them are taken into account the average 
energy dissipation can be as high as 8x 10 6 ergs cm~ 2 s _1 , 
while if one includes areas with lower field strengths 
the average dissipation rate declines, but remains above 
2 x 10 6 ergs cm~ 2 s _1 for the field strengths shown in Fig. 
1101 These values of the average heating are well within 
the observational limits. The total energy dissipated if 
the chromosphere and transition region are also included 
is substantially larger than these numbers, because most 
of the energy is dissipated in the high field strength en- 
vironment of the lower atmosphere. 

3.3. Stratification 
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Fig. 10. — Average heating rate in UV to X-ray emitting gas as Fig. 11. — Surface showing the height of the transition region, 

a function of the underlying photospheric magnetic field strength defined as where the temperature rises above 10 5 K 

threshold 



As limb observation s by TRACE have indi- 
cated for some time ijSchriiver fc McMnllenl HoOfi! 
iTitle fc Schriiverl I2003L private communication), the 
transition region does not occur at a well defined height. 
The variation in transition region height, defined in our 
model as the height where the temperature rises above 
10 5 K, is shown in Fig. for one of the snapshots. 
The height above the photosphere varies from 2.7 Mm 
to 12.3 Mm, with an average of 5.0 Mm, and changes 
over small distances. The intermittency of the height 
of the transition region is a result of the intermittent 
heating, which causes evaporation and re-condensation 
of chromospheric material. 

The stratification of the corona is intermittent on 
even the smallest scales. The temperature is generally 
around 1 MK but there are large deviations from this 
average temperature. There are differences in temper- 
ature of up to 0.7 MK over distances comparable to 
TRACE'S pixel size (~ 375 km). The temperature gra- 
dients are most likely limited by the numerical resolu- 
tion and are therefore only lower limits. Because of 
the large temperature gradients, modelling coronal loops 
as monolithic structures cannot be a good approxima- 
tion J _tlrus_jnakjn^_jnulti^t^ models of loops called 
for IjAschwanden et al.ll2000b|) . 

Figure^Jshows that even at large heights in the corona 
there is gas at low temperatures. Evidence for cool gas 
suspended in the corona has been found through limb 
observations with TRACE based on opacity esti mat es 
(jSchriiver fc McMullenlll999t ISchriiver et alJll999h . The 
cool dense material in the lower part of the corona are 
primarily surges of chromospheric material that flow into 
the corona, or condensations of coronal material on its 
way down. Higher up in the co rona it is typica lly gas un- 
dergoing catastrophic cooling l)Schriiverl l2001'). At times 
surges of material happen simultaneously in both foot 
points of a loop, which brings large amounts of material 
to the top of the loop. If the loop is inclined relative 
to the vertical, accumulated mass at its top can force 
the loop top down and make it almost horizontal along 
part of its length, thus suspending the cold material in 
the corona. An example is the large structure at roughly 
(x, y) = (40, 25) in Fig. 1111 which has a horizontal col- 



umn density of up to 6 x 10 20 cm" 2 at a height of 8.0 Mm, 
making it a visible absorbing structure for UV at the 
limb. The structure is long lived, having existed for 10 
minutes at the end of the simulation sequence. There are 
other similar structures, but too few to give a basis for 
firm estimates of typical characteristics. 

The density has values that vary by at least 2 orders 
of magnitude in the corona, over any horizontal slice. 
The average density and temperature in the corona de- 
pend primarily on the amount of Poynting flux injected 
through the lower boundary and on the choice of the 
chromospheric temperature structure, while the horizon- 
tal variations of density and temperature primarily re- 
flect the variations of properties from loop to loop. In 
general the highest densities are in loop structures with 
low temperatures, persisting all the way to the top of 
the atmosphere. A 2D histogram of the Probability Den- 
sity Function (PDF) of the density at each height for the 
FAL-C chromosphere is displayed in Fig. El which shows 
that the horizontally averaged density stays almost con- 
stant through the part of the corona that we model, with 
an average number density close to 10 9 cm -3 , but with 
values ranging from 10 8 to a few times 10 10 cm -3 , typi- 
cal of or slightly highe r than observed for active regions 
loops (see for instance iMason et alJl!999|) . 

The velocities in the corona are not completely field 
aligned. In the lower corona 80 % of the velocity ampli- 
tude is along the magnetic field, whereas higher in the 
corona the fraction drops to roughly 50 %. In the lower 
corona the strong field regions have velocities that are 
almost fully field aligned, while the lower field strength 
regions have less prevalence for velocities along the mag- 
netic field. In the upper corona the velocities are no 
longer strongly correlated with the magnetic field, and 
the velocities perpendicular to the magnetic field are now 
as strong as those aligned with the magnetic field. In 
general the velocity amplitudes in the corona are high- 
est just above the transition region, where the impulsive 
energy releases can accelerate plasma to high velocities. 
Velocities in this region reach as high as 400 km s _1 but 
with an average of only a few tens of km s . 

3.4. TRACE emission measure 
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Fig. 12. — Histogram of the temperature PDF for each height for 
the FAL-C chromosphere. Dark is high PDF, with a logarithmic 
grey scale. The horizontal average is over plotted (dashed) 



Direct comparison of actual photon count rates be- 
tween TRACE observations and these simulations is 
compl icated for several reasons. iDel Zanna fc Masonl 
( 2003) pointed out that the temperature response func- 
tions of the T RACE filters are based on the Chianti 
v2.0 database ijLandi et alJ I199£^) . which docs not in- 
clude recombinations to and transitions for Fe VIII, and 
this affects the response functions somewhat. Several 
ionisat ion equilibria models are published (see for in- 
stance |Mazzotta et al .l|1998t lArnaud fc Raymond! 11993: 
lArnaud fc Rothenflud Il985|) . with differences that can 
give ch anges in the response funct ion by up to a factor 
of two ijDel Zanna fc Masonll2003l private communica- 
tion). In this wor k, we have used the TRACE response 
function found bv lDel Z anna fc Masonl 112 003). with the 
ionisation states of iMazzo tta et al\ l| 1998D and the ele- 
mental abundances by Fcldman (1992'), calculated at a 
constant electron pressure of 10 15 Kern" 2 . We have also 
made calculations using the local electron pressure in- 
stead of a constant electron pressure, and found that the 
difference is minor (typically less than 10 %), but the 
calculations are much more computationally expensive. 

The images in Fig. El have been produced by folding 
the temperatures and densities of the simulation with the 
response function. In the run with the FAL-C chromo- 
sphere the calculated intensities in the TRACE 171 and 
195 filters are between 60-70 and 15-25 DNs -1 pixel -1 . 
For individual loops chosen from the emulated TRACE 
171 image, that is too high by almost a factor ten for the 
171 filter and a factor 3 for the 195 filter com pared to the 
values given bv IDel Zanna fc Mason (200^). There are 
several possible explanations for this discrepancy. The 
images show a forest of loops of all heights. Some of 
the loops appear broad and diffuse but this is often an 
effect of the perspective — most of these broad loops are 
not connected along the pr ojected axis. The numbers in 
IDel Zanna fc Masonl |2?)03) are for a single, well isolated 
high reaching loop while the large number of loops in 
the emulated images often are hard to isolate. The max- 
imum coronal loop height in the current model is limited 
to ~ 30 Mm; loops penetrating the upper boundary are 
not reliable, since the boundary conditions cannot rigor- 
ously model the effects of the loop top on the part of the 
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Fig. 13. — Histogram of the density PDF for each height for the 
FAL— C chromosphere. Dark is high PDF, with a logarithmic grey 
scale. The horizontal average for the distribution is over plotted 
(dashed) 



loop included in the box. The count rates are therefore 
hard to compare directly. The difference in the count 
rates in the two filters indicates that there is lack of high 
temperature material. Based on the previous work by 
Galsgaard & Nordlund ( 1996), we expect that higher res- 
olution will create more intermittent structures, creating 
a broader range of temperatures and effectively shifting 
gas out of the 171 and 195 filters. 

Cross sections and other proper ties for individual loops 
are ex amined in a separate paper ijGudiksen fc Nordlundl 
12001 . 

4. CONCLUSION 

The increase in computer performance and system 
memory in recent years has made it possible to ap- 
proach the solar coronal heating problem from a new 
angle. Previously mainly qualitative models have been 
developed, with assumptions of the magnetic field struc- 
ture and/or the atmospheric stratification that neces- 
sarily have been simplified. This has kept conclusions 
about the underlying processes indecisive. In recent 
years observational results coming from the SOHO and 
TRACE satellites have been able to put constraints on 
the overall heating function, favouring a heating function 
decreasing exponentially with height l)Schriiver et alJ 
119991 lAschwanden et alJ l2000al 1200 ID . Observations 
of nano- and micro-flares seem to indicate that even 
if one extrapolates the distributions down to energies 
which arc unobserv able, they are not able to supply all 
the energy needed iAs chwandcn & CharbonngajJ l2002t 
lAschwanden fc Parnel]|l2002HParnell fc Jupdl200q) . 

By performing direct numerical simulations we have 
been able to show that moving foot points of the mag- 
netic field around in a way consistent with the observed 
solar photospheric velocity fields inevitably leads to an 
amount of energy dissipation in the corona that is com- 
fortably within the observational limits. The approach 
is nearly ab initio, using only observed facts such as 
the average velocity field properties, a realistic active 
region photospheric magnetic field, a realistic optically 
thin cooling function, Spit zer c onduc tivity, etc. 

As anticipated by iParkerl 1^)72), and as demon- 
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Fig. 14. — Emulated TRACE 171 (left top) and 195 (right top) images both spanning 1.5 box lengths, and the underlying photospheric 
magnetic field 



stra ted numerical l y by iGalseaard Ez Nordlundl l|1996fl 
and lHendrix et all l|1996l ). the magnetic dissipation is ex- 
pected to depend only weakly on the numerical resolu- 
tion. If anything it increases slightly as smaller scales are 
resolved. An aspect that will depend on numerical resolu- 
tion, however, is the time variabil ity and intermittency of 
the he ating. As demonstrated bv lGalsgaard &: Nordlundl 
( 1996), the hierarchy of current sheets in which the mag- 
netic dissipation takes place extends to ever smaller size 
as the numerical resolution is increased. It is presum- 
ably the intrinsic time variability from the reconnection 
events in this hierarchy that gives rise to the flare event 
size power law distribution. From this point of view, 
there is no conflict between the current model and the 
observations that have inspired the nano-flare coronal 
heating models. 

The main results of the coronal modelling is that the 
total dissipated energy is within the observational lim- 
its, making this heating mechanism at least a major 
constituent (and an unavoidable one!) in heating the 
solar corona. The energy release should at least par- 
tially be observable as a dis tribution of flare-like events 
ijGalseaard fc Nordlundl |j"996|) . but we cannot yet esti- 
mate the fraction of the energy released in the form of 
nano-flares, or predict their energy distribution, since 
this requires very high numerical resolution, outside the 
range of present day computer capabilities. 

The magnetic field is in two distinctly different states 
above and below the height where the (3 = 1 layer is 
located. Below, the field is controlled by the photospheric 
gas motions, and is in general in a non-relaxed state. The 
magnetic field there dissipates energy at a high rate, since 
the magnetic energy density is high, the magnetic field 
lines are short, and since oppositely directed field lines 



may be forced towards each other without the magnetic 
pressure having a significant effect on the total pressure. 

In the present simulations roughly 90 % of the to- 
tal dissipated energy is dissipated below the transi- 
tion region. Above the transition region the mag- 
netic field assumes a non-linear near-force-free configu- 
ration, where the dissipated energy is roughly propor- 
tional to the magnetic energy density. Several obser- 
vations point to wards a heating function with precisely 
these properties (IMandrini et alJfeOOOHFolev et al.l | 2002 : 
i Schriiver k, Aschwandenl l2002t [Demoulin et alJ 12003 : 
ISchmieder et alJl2003j) 

The results for the distribution of both temperature 
and density are also well within the observed values, even 
though these values might change slightly depending on 
the stratification of the upper chromosphere. The tem- 
perature is intermittent on the smallest scales, changing 
by as much as 0.7 MK on the scale of a TRACE pixel 
size. The temperature has values from 10 to 3 x 10 6 
K in most of the corona, with a transition zone defined 
as being at temperatures near 5 x 10 5 K varying by as 
much as 9 Mm in height, simila r to what has been in- 
ferred from TRACE observati ons ijSchriiver fc McMuller] 
1999; iTitle &; Schriiverll2003L private communication). 
The average density is almost constant with height in 
the corona, but has values at each height spanning two 
orders of magnitude. 

From the point of view of the energy balance, the den- 
sity and temperature in the corona play different but 
tightly coupled roles. On the one hand the radiative 
cooling function changes relatively little with temper- 
ature, from a few times 10 5 K to a few times 10 6 K, 
while the cooling per unit volume is proportional to the 
square of the density. On the other hand heat conduc- 
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tion scales as T 7 / 2 , but any increase in heat conduction 
also results in more evaporation from the chromosphere, 
so tends to increase the density in the corona. So in this 
sense the temperature is not by itself a sensitive probe for 
the heating in the corona, but has to be taken together 
with density information, and with information about 
the temperature gradients along the magnetic field. 

Making a direct comparison with TRACE observations 
of an active region is hampered by both the question of 
the statistical significance of a single observation, and 
by the insufficient knowledge of the ionisation balance in 
the corona which makes computed photon count rates 
uncertain by up to a factor of two. Furthermore, the 
density in the corona is a function of the stratification in 
the chromosphere, as commented upon earlier. 

In spite of these difficulties the emission measures in 
the TRACE 171 and 195 A bands estimated from the 
simulations are not far from observed values. The dis- 
crepancies (mainly too much emission in the 171 band) 
is most likely due to missing intermittency, in turn due 
to the limited numerical resolution. As discussed in the 
previous Section, increased numerical resolution is ex- 
pected to give a broader distribution of temperatures, 
reaching higher peak values. This will reduce the emis- 
sion measure in the 171 band and will also affect the 195 
band. 

To reproduce the hotter Yohkoh-loops will certainly re- 
quire a temperature distribution that extends to higher 
values than what we find here. Increased heating will in- 
evitably result if a more complete representation of the 
solar velocity field is used. In the current simulations 
there is a resolution cut-off in the velocity driving to- 
wards smaller scales, and a box size cut-off towards larger 
scales. It is interesting to note that scaling expressions 
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